Explore genetic evolution only:

Constant Environment

To explore this we used the following parameters space :

constant_env=getAlllSummaries("unifPop50k")

I was not keeping track of the outpute rate so the convertion in real generation has to be done by hand:

constant_env$outputrate=10
constant_env=updateScale(constant_env)

Trajectories

Check the trajectories of the simulation for the population size for different mutation rate \(\mu\)(from left to right: \(\mu = 0, \mu = 0.001,\mu=0.01\)), different selective pressure \(\sigma_s\) (from top to bottomright: \(\sigma_s = 0.1, \sigma_smu = 0.2,\sigma_s=0.4,\sigma_s=1000\))and different carrying capacity \(K\)

#pdf("trajectory_N_const.pdf",width=12,height=16)
plotAllTrajVar(constant_env,m=0.05,E=0,obs="N")

#dev.off()
#pdf("trajectory_var_x_const.pdf",width=12,height=16)
plotAllTrajVar(constant_env,m=0.05,E=0,obs="var_x",ylim=c(0,.005))

#dev.off()

Variance at equilibrium

Check the variance at the equilibrium:

#pdf("equilibrium_var_x_const.pdf",width=12,height=16)
plotAllVariableSummaries(constant_env,E=0,estimate=eq2833b) #we use only simulations without noise

#dev.off()

Constant but noisy environment

we use different \(delta\) to generate environment with noise of growing variance

Environment

d=c(0.00,0.01,0.10,0.20,0.40,1.00,10.00)
cols=colorRampPalette(c("black","grey"))(length(d))
names(cols)=d
omegas=1:3

plot(1,1,xlim=c(1,500),ylim=c(-10,10),type="n",xlab="t",ylab=expression(theta[t]))
     for ( i in rev(d))
         lines(1:500,environment(500,omega=0,delta=i),col=cols[as.character(i)])
legend("topleft",legend=paste0("v=",d),col=cols,lty=1,bg="white")

To explore this we used the following parameters space:

v=getAlllSummaries("growingDelta",exclude="outputrate")
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
noisy_env = rbind(v,constant_env[,colnames(v)])
noisy_env = noisy_env[ noisy_env$mu >0,]

Trajectories

Check the trajectory off the effective population size:

#pdf("trajectory_N_noise.pdf",width=12,height=16)
plotAllTrajVar(noisy_env,m=0.05,E=0,obs="N")

#dev.off()

Check the trajectory off the mean variance size:

#pdf("trajectory_var_x_noise.pdf",width=12,height=16)
plotAllTrajVar(noisy_env,m=0.05,E=0,obs="var_x",ylim=c(0,.005))

#dev.off()

Variance at equilibrium

Check and compare the final variance :

#pdf("equilibrium_var_x_noise.pdf",width=12,height=16)
plotAllVariableSummaries(noisy_env,E=0,estimate=eq2833b) 

#dev.off()

Some of the data seems to be missing because when \(\delta>\sigma\) we have extinctions, has shown in the graph:

Extinctions

#pdf
sigmas=unique(noisy_env$sigma)
cols=colorRampPalette(c("darkgreen","yellow"))(length(sigmas))
names(cols)=sigmas
tmpc=noisy_env[ noisy_env$sigma==0.4 & noisy_env$mu == 0.01 & noisy_env$m == 0.2 ,]
plot(tmpc$delta,tmpc$extinction,log="xy",type="n",ylab="extinction time",xlab=expression(delta),main=bquote(mu == 0.01 ~ omega == 0))
for(s in sigmas){
    tmpc=noisy_env[ noisy_env$sigma==s & noisy_env$mu == 0.01 & noisy_env$m == 0.2 ,]
    Ks=sort(unique(tmpc$K))
    for(k in 1:length(Ks)){
        t=tmpc[tmpc$K==Ks[k],]
        r=tapply(t$extinction,t$delta,mean)
        lines(sort(unique(t$delta)),r,type="b",pch=k,col=cols[as.character(s)])
    }
}
        legend("right",
               legend=c(paste0("K=",Ks),sapply(sigmas,function(s)as.expression(bquote(sigma==.(s))))),
               col=c(rep(1,length(Ks)),cols),
               lty=c(rep(NA,length(Ks)),rep(1,length(sigmas))),
               pch=c(seq_along(Ks),rep(NA,length(sigmas)))
               )

#dev.off()

Constant, noisy, with autocorrelation

Environment

We introduce autocorrelation in the environment. To introduce autocorrelation we increase slightly \(omega\)

par(mfrow=c(1,3))
plot(1:500,environment(500,omega=1,delta=.1),type="l",xlab="t",ylab=expression(theta[t]))
plot(1:500,environment(500,omega=2,delta=.1),type="l",xlab="t",ylab=expression(theta[t]))
plot(1:500,environment(500,omega=3,delta=.1),type="l",xlab="t",ylab=expression(theta[t]))

For this we re-run previous experiment but with \(\omega \in {1,2}\)

omega1=getAlllSummaries(c("omega1growingDelta"))
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
omega1=updateScale(omega1)
omega1$omega=1

omega2=getAlllSummaries(c("omega2growingDelta"))
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
omega2=updateScale(omega2)
omega2$omega=2

noisy_env_omega = rbind(v,constant_env[,colnames(v)],omega1[,colnames(v)],omega2[,colnames(v)])

Extinctions

Let’s look at the extinction again for the different omegas

par(mfrow=c(1,3))
omegas=unique(noisy_env_omega$omega)
sigmas=sort(unique(noisy_env_omega$sigma))
deltas=sort(unique(noisy_env_omega$delta))
Ks=sort(unique(noisy_env_omega$K))
for(o in omegas){
    cols=colorRampPalette(c("darkgreen","yellow"))(length(sigmas))
    names(cols)=sigmas
    tmpc=noisy_env_omega[ noisy_env_omega$sigma==0.4 & noisy_env_omega$mu == 0.01 & noisy_env_omega$omega == o & noisy_env_omega$m == 0.2 ,]
    xrange=range(noisy_env_omega$delta[noisy_env_omega$omega>0])
    plot(tmpc$delta,tmpc$extinction,log="xy",type="n",ylab="extinction time",xlab=expression(delta),main=bquote(mu == 0.01 ~ omega == .(o)),xlim=xrange)
    for(s in sigmas){
        tmpc=noisy_env_omega[ noisy_env_omega$sigma==s & noisy_env_omega$mu == 0.01 & noisy_env_omega$m == 0.2 & noisy_env_omega$omega == o,]
        for(k in 1:length(Ks)){
            t=tmpc[tmpc$K==Ks[k],]
            r=tapply(t$extinction,t$delta,mean)
            lines(sort(unique(t$delta)),r,type="b",pch=k,col=cols[as.character(s)])
        }
    }
    legend("bottomleft",
           legend=c(paste0("K=",Ks),sapply(sigmas,function(s)as.expression(bquote(sigma==.(s))))),
           col=c(rep(1,length(Ks)),cols),
           lty=c(rep(NA,length(Ks)),rep(1,length(sigmas))),
           pch=c(seq_along(Ks),rep(NA,length(sigmas)))
           )
}

Trajectories

We can thus reproduce he trajectories and variance equilibrium using simulation with \(\omega=1\) et \(\omega=2\)

Check the trajectory off the effective population size:

#pdf("trajectory_N_omega1.pdf",width=12,height=16)
plotAllTrajVar(omega1,m=0.05,E=0,obs="N")

#dev.off()
#pdf("trajectory_N_omega2.pdf",width=12,height=16)
plotAllTrajVar(omega2,m=0.05,E=0,obs="N")

#dev.off()

Check the trajectory off the mean variance size:

#pdf("trajectory_var_x_omega1.pdf",width=12,height=16)
plotAllTrajVar(omega1,m=0.05,E=0,obs="var_x",ylim=c(0,.002))

#dev.off()
#pdf("trajectory_var_x_omega2.pdf",width=12,height=16)
plotAllTrajVar(omega2,m=0.05,E=0,obs="var_x",ylim=c(0,.002))

#dev.off()

Variance at equilibrium

And check the variance at the end of the equilibrium:

constant_env$omega=0
omega1_wdz=rbind(omega1,constant_env[constant_env$mu>0,])
omega2_wdz=rbind(omega2,constant_env[constant_env$mu>0,])
omega1_wdz=omega1_wdz[omega1_wdz$sigma<10,]
omega2_wdz=omega2_wdz[omega2_wdz$sigma<10,]
#pdf("equilibrium_var_x_omega1.pdf",width=12,height=16)
plotAllVariableSummaries(omega1_wdz,E=0,estimate=eq2833b) #we use only simulations without noise

#dev.off()
#pdf("equilibrium_var_x_omega2.pdf",width=12,height=16)
plotAllVariableSummaries(omega2_wdz,E=0,estimate=eq2833b) #we use only simulations without noise

#dev.off()

Moving optimum \(\theta_t=vt\)

Environment

par(mfrow=c(3,3))
deltas=c(0,0.1,1)
vts=c(0.001,0.002,0.003)
cols=colorRampPalette(c("darkgreen","green"))(length(vts))
names(cols)=vts
omegas=1:3

for(d  in deltas){
for(o  in omegas){
plot(1,1,xlim=c(1,500),ylim=c(-(d+d*.1),1.1*500*max(vts)+d),type="n",xlab="t",ylab=expression(theta[t]),main=bquote(delta == .(d)~omega == .(o)))
     for ( vt in vts)
         lines(1:500,environment(500,omega=o,delta=d,vt=vt),col=cols[as.character(vt)])
legend("topleft",legend=paste0("v=",vts),col=cols,lty=1)
}
}

moving_theta=getAlllSummaries("movingTheta")
moving_theta=updateScale(moving_theta)

Extinctions

#pdf("extinctionsMovingTheta.pdf",width=17,height=4)
par(mfrow=c(1,3))
rates=sort(unique(moving_theta$vt))
sigmas=sort(unique(moving_theta$sigma))
deltas=sort(unique(moving_theta$delta))
Ks=sort(unique(moving_theta$K))
cols=colorRampPalette(c("darkgreen","yellow"))(length(sigmas))
    names(cols)=sigmas
for(d in deltas){
    tmpcs=moving_theta[  moving_theta$mu == 0.001 & moving_theta$delta == d & moving_theta$m == 0.1 ,]
    #xrange=range(moving_theta$delta)
    plot(tmpcs$vt,tmpcs$extinction,log="y",type="n",ylab="extinction time",xlab=expression(v),main=bquote(delta == .(d)~mu == 0.01 ~  omega == .(o)))
    for(s in sigmas){
        tmpc=tmpcs[ tmpcs$sigma==s ,]
        for(k in 1:length(Ks)){
            t=tmpc[tmpc$K==Ks[k],]
            r=tapply(t$extinction,t$vt,mean)
            lines(sort(unique(t$vt)),r,type="b",pch=k,col=cols[as.character(s)])
        }
    }
    legend("bottomleft",
           legend=c(paste0("K=",Ks),sapply(sigmas,function(s)as.expression(bquote(sigma==.(s))))),
           col=c(rep(1,length(Ks)),cols),
           lty=c(rep(NA,length(Ks)),rep(1,length(sigmas))),
           pch=c(seq_along(Ks),rep(NA,length(sigmas)))
           )
}

#dev.off()

Trajectories

Check the trajectory of the effective population size:

for(t in seq_along(rates)){
#pdf(paste0("trajectory_N_moving",t,".pdf"),width=12,height=16)
plotAllTrajVar(moving_theta[moving_theta$vt == rates[t],],m=0.1,E=0,obs="N")
#dev.off()
}

Check the trajectory of the mean variance size:

for(t in seq_along(rates)){
#pdf(paste0("trajectory_var_x_moving",t,".pdf"),width=12,height=16)
plotAllTrajVar(moving_theta[moving_theta$vt == rates[t],],m=0.1,E=0,obs="var_x",ylim=c(0,.005))
#dev.off()
}

Check the trajectory of the distance to optimum:

for(t in seq_along(rates)){
#pdf(paste0("trajectory_dist_moving",t,".pdf"),width=12,height=16)
plotAllTrajVar(moving_theta[moving_theta$vt == rates[t],],m=0.1,E=0,obs="dist")#,ylim=c(0,.005))
#dev.off()
}

Variance at equilibrium

Check and compare the final variance :

for(t in seq_along(rates)){
#pdf(paste0("equilibrium_var_x_moving",t,".pdf"),width=12,height=16)
plotAllVariableSummaries(moving_theta[moving_theta$vt == rates[t],],E=0,estimate=eq2833b) 
#dev.off()
}